namespace Dune {
  namespace Compose {
    template<typename T, typename X>
    class FastEigenvalueDGF : public Dune::PDELab::GridFunctionBase<
      Dune::PDELab::GridFunctionTraits<
        typename T::Traits::GridViewType,
        typename T::template Child<0>::Type::Traits::FiniteElementType
          ::Traits::LocalBasisType::Traits::RangeFieldType,
        T::template Child<0>::Type::Traits::FiniteElementType
          ::Traits::LocalBasisType::Traits::dimRange,
        typename T::template Child<0>::Type::Traits::FiniteElementType
          ::Traits::LocalBasisType::Traits::RangeType>,
      FastEigenvalueDGF<T,X> >
    {
      typedef T GFS;

      typedef Dune::PDELab::GridFunctionBase<
        Dune::PDELab::GridFunctionTraits<
          typename T::Traits::GridViewType,
          typename T::template Child<0>::Type::Traits::FiniteElementType
            ::Traits::LocalBasisType::Traits::RangeFieldType,
          T::template Child<0>::Type::Traits::FiniteElementType
            ::Traits::LocalBasisType::Traits::dimRange,
          typename T::template Child<0>::Type::Traits::FiniteElementType
            ::Traits::LocalBasisType::Traits::RangeType>,
        FastEigenvalueDGF<T,X> > BaseT;

     typedef Dune::FieldVector<
        typename T::template Child<0>::Type::Traits::FiniteElementType
          ::Traits::LocalBasisType::Traits::RangeFieldType,
        TypeTree::StaticDegree<T>::value
        > FV;
     typedef typename T::template Child<0>::Type::Traits::FiniteElementType
       ::Traits::LocalBasisType::Traits::RangeType RT;

    public :
      typedef typename BaseT::Traits Traits;
      typedef typename T::template Child<0>::Type ChildType;
      typedef typename ChildType::Traits::FiniteElementType
        ::Traits::LocalBasisType::Traits LBTraits;
      typedef typename LBTraits::RangeFieldType RF;

      template<typename EQ>
      FastEigenvalueDGF(const GFS& gfs, const X& x_, EQ& eq)
        : pgfs(stackobject_to_shared_ptr(gfs))
        , lfs(gfs)
        , lfs_cache(lfs)
        , x_view(x_)
        , xl(gfs.maxLocalSize())
      {
        this->eq = &eq;
      }

      template<typename EQ>
      FastEigenvalueDGF (std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x_, EQ& eq)
        : pgfs(gfs)
        , lfs(*gfs)
        , lfs_cache(lfs)
        , x_view(*x_)
        , xl(pgfs->maxLocalSize())
      {
        this->eq = &eq;
      } 

      template <typename EQ> 
      void set_equation (EQ& eq) {
        this->eq = &eq;
      }

      inline void evaluate(const typename Traits::ElementType& e,
                           const typename Traits::DomainType& x,
                           typename Traits::RangeType& y) const
      {
        // get and bind local functions space
        lfs.bind(e);
        lfs_cache.update();
        x_view.bind(lfs_cache);
        x_view.read(xl);
        x_view.unbind();
        y = 0.0;
        FV yv;
        double yvtmp[TypeTree::StaticDegree<T>::value];
        // std::cout << "in getdt, show dimRange: " << T::template Child<0>::Type::Traits::FiniteElementType
        //   ::Traits::LocalBasisType::Traits::dimRange << std::endl;

        // std::cout << "in getdt, show staticdgree: " << TypeTree::StaticDegree<T>::value << std::endl;
        std::vector<RT> yb;
        // double yv[10];
        // loop over VectorGFS and calculate k-th derivative of k-th child
        for(unsigned int k=0; k != TypeTree::degree(lfs); ++k) {
        // for(unsigned int k=0; k != 1; ++k) {
          yv[k] = 0.0;
          lfs.child(k).finiteElement().localBasis().evaluateFunction(x,yb);
          for(typename LFS::Traits::SizeType i=0; i<lfs.child(k).size(); i++) {
              yv[k] += xl[lfs.child(k).localIndex(i)]*yb[i];
          }
          yvtmp[k] = yv[k];
          // std::cout <<  "yv["<< k <<"]: " << yv[k] <<std::endl;
          // y = yv[0]+yv[1]+yv[2]+yv[3];
        }
        //y = yv[0]; //ok
        // y = 108; //ok;
        // y = yvtmp[0]+yvtmp[1]+yvtmp[2]+yvtmp[3];
        // y = yv[0]+yv[1]+yv[2]+yv[3];
        y = eq->local_max_eigenvalue(yvtmp);
      }

      inline const typename Traits::GridViewType& getGridView() const {
        return pgfs->getGridView();
      }

    private :
      typedef Dune::PDELab::LocalFunctionSpace<GFS> LFS;
      typedef Dune::PDELab::LFSIndexCache<LFS> LFSCache;
      typedef typename X::template ConstLocalView<LFSCache> XView;
      std::shared_ptr<GFS const> pgfs;
      mutable LFS lfs;
      mutable LFSCache lfs_cache;
      mutable XView x_view;
      mutable std::vector<RF> xl;
      std::shared_ptr<const X> px;
      Dune::Equation::EquationBase* eq;
    };

//    template<typename RF>
//    void Suggestdt(RF& maxfastev) 
//    {
//      maxfastev = 1.0;
//    };

    template<typename GV, typename GF, typename RF>
    void maxfasteigenvalue(const GV& gv, const GF& gf, RF& ymax)
    {
      typedef typename GF::Traits::DomainType DT;
      typedef typename GF::Traits::RangeType RT;
      DT x(0.0);
      RT y;
      typedef typename GV::Traits::template Codim<0>::Iterator ElementIterator;
      // FIXME: at the end of simulation, this loop lead to uncorrect stop behaviour when use mult-nodes.
      for (ElementIterator it = gv.template begin<0>();
        it!=gv.template end<0>(); ++it) {
        gf.evaluate(*it,x,y);
        if (ymax <= std::abs(y)) {
          ymax = std::abs(y);
        }
      }
      // return ymax;
    };

    template<typename RF>
    void getdt(RF& maxfastev, RF& dt,Dune::ParameterTree& ptree) {
      // cfl condition: max_fast_eigenvalue dt/dx
      dt = ptree.get<double>("fem.cfl") * 
        (ptree.get<double>("grid.structured.LX") / 
        ptree.get<int>("grid.structured.NX")) / maxfastev;
    };
  } // Compose
} // Dune
